Model selection for cluster data analysis

ABSTRACT

A model selection method is provided for choosing the number of clusters, or more generally the parameters of a clustering algorithm. The algorithm is based on comparing the similarity between pairs of clustering runs on sub-samples or other perturbations of the data. High pairwise similarities show that the clustering represents a stable pattern in the data. The method is applicable to any clustering algorithm, and can also detect lack of structure. We show results on artificial and real data using a hierarchical clustering algorithm.

FIELD OF THE INVENTION

The present invention relates to the use of learning machines toidentify relevant patterns in datasets containing large quantities ofdiverse data, and more particularly to a method and system forunsupervised learning for determining an optimal number of data clustersinto which data can be divided to best enable identification of relevantpatterns.

BACKGROUND OF THE INVENTION

Knowledge discovery is the most desirable end product of datacollection. Recent advancements in database technology have lead to anexplosive growth in systems and methods for generating, collecting andstoring vast amounts of data. While database technology enablesefficient collection and storage of large data sets, the challenge offacilitating human comprehension of the information in this data isgrowing ever more difficult With many existing techniques the problemhas become unapproachable. Thus, there remains a need for a newgeneration of automated knowledge discovery tools.

As a specific example, the Human Genome Project has completed sequencingof the human genome. The complete sequence contains a staggering amountof data, with approximately 31,500 genes in the whole genome. The amountof data relevant to the genome must then be multiplied when consideringcomparative and other analyses that are needed in order to make use ofthe sequence data. As an illustration, human chromosome 20 alonecomprises nearly 60 million base pairs. Several disease-causing geneshave been mapped to chromosome 20 including various autoimmune diseases,certain neurological diseases, type 2 diabetes, several forms of cancer,and more, such that considerable information can be associated with thissequence alone.

One of the more recent advances in determining the functioningparameters of biological systems is the analysis of correlation ofgenomic information with protein functioning to elucidate therelationship between gene expression, protein function and interaction,and disease states or progression. Proteomics is the study of the groupof proteins encoded and regulated by a genome. Genomic activation orexpression does not always mean direct changes in protein productionlevels or activity. Alternative processing of mRNA or posttranscriptional or post-translational regulatory mechanisms may causethe activity of one gene to result in multiple proteins, all of whichare slightly different with different migration patterns and biologicalactivities. The human proteome is believed to be 50 to 100 times largerthan the human genome. Currently, there are no methods, systems ordevices for adequately analyzing the data generated by such biologicalinvestigations into the genome and proteome.

Clustering is a widely used approach for exploratory data analysis.Clustering analysis is unsupervised learning—it is done withoutsuggestion from an external supervisor; classes and training examplesare not given a priori. The objective of clustering is to group datapoints into “meaningful” subsets. There is no agreed upon definition ofthe clustering problem, and various definitions appear in theliterature. For example, clustering has been defined as a search forsome “natural” or “inherent” grouping of the data. However, mostclustering algorithms do not address this problem. The vast majority ofclustering algorithms produce as their output either a dendogram or apartition into a number of clusters, where the number of clusters iseither the input, or there is some other parameter(s) that controls thenumber of clusters. In either case, a model selection technique isrequired in order to choose the model parameter, or in the case ofhierarchical algorithms, to determine which level of the dendogramrepresents the “inherent” structure of the data.

A few examples of applications of clustering include (1) analysis ofmicroarray data, where co-expressed genes are found, and the assumptionis that co-expression might be a sign of co-regulation; (2) in medicaldatasets (gene expression data, clinical data etc.), where patients aredivided into categories; (3) in any set or set of measurements to detecttrends or artifacts in the measurement protocol; and (4) in informationretrieval to partition text according to categories.

Most clustering algorithms either produce a hierarchical partitioning ofthe data into smaller and smaller clusters, or produces a partition of adataset into a number of clusters that depend on some input parameter(the number of clusters or some other parameter(s)). The questionremains, however, of how to set the input parameter, or how to determinewhich level of the tree representation of the data to look at:Clustering algorithms are unsophisticated in that they provide noinsight into the level of granularity at which the “meaningful” clustersmight be found. Occasionally, there may be prior knowledge about thedomain that facilitates making such a choice. However, even in suchcases, a method for determining the granularity at which to look at thedata is required. This is seen as the problem of finding the optimalnumber of clusters in the data, relative to some clustering algorithm.

E. Levine and E. Domany in “Resampling Method for UnsupervisedEstimation of Cluster Validity”, Neural Comp. 13, 2573-2593 (2001),assign a figure of merit to a clustering solution according to itssimilarity to clusterings of subsamples of the data. The “temperature”parameter of their clustering algorithm is selected according to amaximum of the similarity measure. However, in real data, such a maximumdoes not often occur. Other model selection techniques have difficultydetecting the absence of structure in the data, i.e., that there is asingle cluster. Further, many of algorithms make assumptions as tocluster shape, and do not perform well on real data, where the clustershape is generally not known. Accordingly, other methods for clusteringare needed. The present invention is directed to such a method.

SUMMARY OF THE INVENTION

In an exemplary embodiment, a system and method are provided foranalysis of data by grouping the input data points into meaningfulsubsets of data. The exemplary system comprises a storage device forstoring at least one data set, and a processor for executing theclustering algorithm. In the preferred embodiments, the clusteringalgorithm is a k-means or hierarchical clustering algorithm.

In the method of the present invention, a “good granularity level”,i.e., an optimal number of clusters, is defined as one at which aclustering solution is stable with respect to some perturbation of thedata, such as noise or sub-sampling. For each level of granularity, thealgorithm chooses a number of pairs of sub-samples or otherperturbations, clusters each subsample with the chosen level ofgranularity, then computes a similarity between pairs of clusteringsolutions. For each level of granularity, the distribution of thesimilarity between pairs of clustering solutions is computed, then thehighest level of granularity for which highly similar solutions areobtained under perturbation is chosen.

According to the present invention, the probability distribution of thesimilarity measure is evaluated in order to find the “true” number ofclusters based on a similarity measure between clustering solutions. Adot product between pairs of clustering solutions is normalized toprovide a similarity measure. Other similarity measures proposed in theliterature can also be expressed in terms of this dot product.

In one embodiment, the inventive method provides means for extraction ofinformation from gene expression profiles, in particular, extraction ofthe most relevant clusters of temporal expression profiles from tens ofthousands of measured profiles. In one variation of the presentembodiment, the clustering algorithm is the k-means algorithm. Otherembodiments include all pairwise clustering methods (e.g., hierarchicalclustering) and support vector clustering (i.e., a support vectormachine (SVM)). In one variation, the fit involves affinetransformations such as translation, rotation, and scale. Othertransformations could be included, such as elastic transformations. Inthe case of k-means, the algorithm is applied in a hierarchical mannerby sub-sampling the data and clustering the sub-samples. The resultingcluster centers for all the runs are then clustered again. The resultingcluster centers are considered to be the most significant profiles. Tofacilitate the clustering task, a ranking of the genes is firstperformed according to a given quality criterion combining saliency(significant difference in expression in the profile), smoothness, andreliability (low noise level). Other criteria for ranking include thelocal density of examples.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1 a-1 c are a plot of data which is a mixture of four Gaussians, ahistogram of the correlation score for the Gaussian data mixture, and aplot showing the overlay of cumulative distributions of the correlationscore, respectively.

FIGS. 2 a and 2 b are histogram of the correlation score for the yeastgene expression data and an overlay of the cumulative distributions ofthe correlation score, respectively.

FIGS. 3 a and 3 b are a histogram of the correlation score for 208points uniformly distributed on the unit cube and an overlay of thecumulative distributions of the correlation score, respectively.

FIGS. 4 a-4 c are a plot of the first two principle components, thehistogram of the correlation score, and an overlay of the cumulativedistributions of the correlation score; k=2 (rightmost) to k−10(leftmost), respectively.

FIG. 5 is a graphic representation of gene pre-selection.

FIG. 6 is a graph of the average profiles of the clusters of the 8clustering runs and their grouping into meta-clusters.

FIG. 7 is a graph of curve fitting with affine transformations.

FIG. 8 is a graph of curve fitting with affine transformations,clustering obtained from the top 800 best-quality genes.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

A typical computer system for running the inventive clustering algorithmis a Pentium®-class processor with an interface for inputting thedataset(s) of interest and a memory for storing the dataset. Alternativeprocessors include configurable logic processors made up of SRAM-basedfield programmable gate arrays (FPGAs) or configurable system on chip(CSOC) architectures, which include a processor and an array ofconfigurable logic cells on a single chip, either of which can be usedto accelerate computer-intensive operations such as clusteringalgorithms on large datasets. In addition, parallel implementations inmulti-computers have been used for executing clustering algorithms forextracting knowledge from large-scale data repositories. Selection of anappropriate computer system for executing the inventive clusteringmethod is within the level of skill in the art.

The clustering model selection algorithm works with the help of ascoring function that provides a similarity measure between twolabelings.

Let X={x₁, . . . , x_(n)}, and x_(i)∈R^(d) be the dataset to beclustered. A labeling L is a partition of X into k subsets S₁, . . . ,S_(k). It can be represented by a function c: X→{1, . . . , k} wherec(x_(i)) is the cluster to which x_(i) belongs.

A less compact representation of a labeling, which may be useful is thefollowing representation by a matrix C with components: $\begin{matrix}{C_{ij} = \left\{ \begin{matrix}1 & {{{if}\quad x_{i}\quad{and}\quad x_{j}\quad{belong}\quad{to}\quad{the}\quad{same}\quad{cluster}},} \\0 & {otherwise}\end{matrix} \right.} & (1)\end{matrix}$

Note that this representation is label independent, i.e. there is noneed to assign each point a label in {1, . . . , k}. This may be viewedas an adjacency matrix, where each cluster is a connected component ofthe graph. This representation can also be converted into arepresentation of soft cluster labeling.

Let labelings L₁ and L₂ have matrix representations C⁽¹⁾ and C⁽²⁾respectively, so that $\begin{matrix}{\left( {\mathcal{L}_{1},\mathcal{L}_{2}} \right) = {\sum\limits_{i,j}^{\quad}\quad{C_{ij}^{(1)}C_{ij}^{(2)}}}} & (2)\end{matrix}$

This dot product computes the number of pairs of vectors clusteredtogether, and can also be interpreted as the number of common edges ingraphs represented by C⁽¹⁾ and C⁽²). A naïve method for computing thedot product by going over all pairs of points has complexity O(n²).However, it can be computed in linear time.

As a dot product,

L₁, L₂

satisfies the Cauchy-Schwartz inequality:

L₁, L₂

≦{square root}{square root over (

L₁, L₁

L₂, L₂

)}. The correlation score, which is a normalized version of the dotproduct, is: $\begin{matrix}{{{cor}\left( {\mathcal{L}_{1},\mathcal{L}_{2}} \right)} = \frac{\left\langle {\mathcal{L}_{1},\mathcal{L}_{2}} \right\rangle}{\sqrt{\left\langle {\mathcal{L}_{1},\mathcal{L}_{1}} \right\rangle\left\langle {\mathcal{L}_{2},\mathcal{L}_{2}} \right\rangle}}} & (3)\end{matrix}$

A. K. Jain and R. C. Dubes, Algorithms for clustering data (PrenticeHall, Englewood Cliffs, N.J., 1988) provide a number of scores forcomparing a labeling produced by a clustering algorithm with a “goldstandard” labeling. This technique is known as “external” validation. Incontrast, the present invention proposes the use of scoring functionsfor “internal” validation that do not require the “gold standard”. Inthe following, it is shown that two commonly used scoring functions canbe expressed in terms of the dot product defined above in Equation 2.Given two matrices C⁽¹⁾, C⁽²⁾ with 0-1 entries, let N_(ij) i,j∈{0, 1} bethe number of entries on which C⁽¹⁾ and C⁽²⁾ have values i and j,respectively. Define the matching coefficient as the fraction of entrieson which the two matrices agree: $\begin{matrix}{{M\left( {\mathcal{L}_{1},\mathcal{L}_{2}} \right)} = {\frac{N_{00} + N_{11}}{N_{00} + N_{01} + N_{10} + N_{11}}.}} & (4)\end{matrix}$

The Jaccard coefficient is the corresponding ratio when “negative”matches are ignored: $\begin{matrix}{{J\left( {\mathcal{L}_{1},\mathcal{L}_{2}} \right)} = {\frac{N_{11}}{N_{01} + N_{10} + N_{11}}.}} & (5)\end{matrix}$

The Jaccard coefficient is more appropriate when the clusters arerelatively small, since in that case, the N₀₀ term will be the dominantfactor even if the solution is far from the true one. These scores canbe expressed in terms of the labeling dot product and the associatednorm according to the following proposition:

Let C⁽¹⁾ and C⁽²⁾ be the matrix representations of labelings L₁ and L₂,respectively. Then: $\begin{matrix}{{J\left( {\mathcal{L}_{1},\mathcal{L}_{2}} \right)} = \frac{{< C^{(1)}},{C^{(2)} >}}{{< C^{(1)}},{C^{(1)} > {+ {< C^{(2)}}}},{C^{(2)} > {- {< C^{(1)}}}},{C^{(2)} >}}} \\{{M\left( {\mathcal{L}_{1},\mathcal{L}_{2}} \right)} = {1 - {\frac{1}{n^{2}}{{C^{(1)} - C^{(2)}}}^{2}}}}\end{matrix}$

This is a result of the observation that N₁₁=

C⁽¹⁾, C⁽²⁾

, N₀₁=

1_(n)−C⁽¹⁾,C⁽²⁾

, N₁₀=

C⁽¹⁾, 1_(n)−C⁽²⁾

, N₀₀=

1_(n)−C⁽¹⁾, 1_(n)−C⁽²⁾

, where 1_(n) is an n×n matrix with entries equal to 1. The aboveexpression for the Jaccard coefficient shows that it is close to thecorrelation score.

The following list provides the routine for the model explorer algorithmaccording to the present invention:

-   -   Input: A dataset D, k_(min), k_(max)    -   Require: A clustering algorithm, cluster (D, k)    -   A scoring function for label comparison, s(L₁, L₂)        -   f=0.8        -   for k=k_(min) to k_(max) do            -   for i=1 to maximum iterations do                -   sub₁=subsamp(D,f) {a sub-sample with fraction f of                    the data}                -   sub₂=subsamp(D,J)=                -   L₁=cluster (sub₁, k)                -   L₂=cluster (sub₂, k)                -   Intersect=sub_(1∩)sub₂                -   Score (i, k)=s(L₁(Intersect), L₂(Intersect))                    {Compute the score on the intersection of the two                    labels}            -   end for        -   end for

When one looks at a cloud of data points, and at a subsample of it for asampling ratio, f (fraction of points sampled) is not much smaller than1 (say f>0.5), one usually observes the same general structure. Thus, itis also reasonable to postulate that a clustering algorithm has capturedthe inherent structure in a dataset if clustering solutions overdifferent subsamples are similar, e.g., according to one of thesimilarity measures introduced in the previous section. Thus, “inherentstructure” is defined as structure that is stable under sub-sampling (oralternatively, perturbing the data). Given a clustering algorithm and adataset, this translates into a search for the best number of clustersto use in the particular instance. Note that one can also extend thesearch to look for a set of features where structure is apparent,however, in this case, the same set of features is kept.

The inventive clustering algorithm receives as input a dataset (orsimilarity/dissimilarity matrix) and a parameter k that controls eitherdirectly or indirectly the number of clusters that the algorithmproduces. This convention is applicable to hierarchical clusteringalgorithms as well: given k, the tree is cut so that k clusters areproduced. Next, characterize the stability of the clustering for each k.This is accomplished by producing a set of clusterings of sub-samples ofthe data, and comparing the labels of the intersection of pairs ofclusterings using, for example, the correlation similarity measure. Thisis performed for increasing values of k (see above for details). Thedistribution of the scores for the different values of k is thencompared (see, e.g., FIGS. 1 a-1 e). The idea is that when there isstructure in the data that is well described by the clustering algorithmwith that particular value of k, many sub-samples will produce similarclusterings, and their pairwise similarity score will be concentratedclose to 1.

Each sub-sample contains a fixed fraction of the data, f The actualsubsampling can be implemented in various ways:

1. Select each sample independently so that the size of the intersectionbetween two samples is random.

2. Select together pairs of samples by first selecting theirintersection, then selecting the rest of the data to complete thefraction f.

3. Fix one clustering solution to be one produced on the whole dataset(This third option was used by Levine and Domany (supra) to give afigure of merit to a particular clustering solution.)

For k=1, all clusterings are the same. This also holds fork=n, where nis the number of data points; in this case every point is in a differentcluster. When the number of clusters becomes large so that there is asmall number of points in each cluster, the solution becomes stable. Thevalue off should not be too low so that there not all clusters arerepresented in a sub-sample. In the Examples provided below, the shapeof the distribution did not depend very much on the specific value offAny value between 0.6 and 0.9 worked well.

EXAMPLES

In this section experiments on artificial and real data are described.In all the experiments the distribution of the correlation score isshown. Equivalent results were obtained using other scores as well. Theparameter values f=0.8 and 200 pairs of solutions were compared for eachk. A hierarchical clustering algorithm was used, with the Ward criterionfor merging clusters (see, e.g., Jain and Dubes, supra). Similar resultswere obtained using other hierarchical clustering methods (complete andaverage linkage). The advantage of using hierarchical clustering methodsis that the same set of clusterings can be used for all values of k.

Example 1 Gaussian Data

Referring first to FIGS. 1 a-1 c, FIG. 1 a shows a mixture of fourGaussians. The histograms of the score for varying values of k for thisdata is plotted in FIG. 1 b. Histograms are shown for each value of k inthe range of 2 to 7. Observations regarding the histograms are that atk=2, there is a peak at 1, since almost all the runs discriminatedbetween the two upper and two lower clusters. At k=3, most runsseparated the two lower clusters, and at k=4 most runs found the“correct” clustering as is reflected in the distribution of scores thatis still close to 1.0. At k>4 there is no longer essentially onepreferred solution. There is, in fact, a wide variety of solutions,evidenced by the widening spectrum of the similarities. FIG. 1 c plotsthe cumulative distributions of the correlation score for each k, wherek=2 at the rightmost side of the plot (at peak 1), and k=7 being theleftmost curve.

Example 2 DNA Microarray Data

The next dataset considered was the yeast DNA microarray data of M.Eisen et al. (“Genetics cluster analysis and display of genome-wideexpression patterns”, Proc. Natl. Acad. Sci. USA, 95: 14863-14868,December 1998. The data is a matrix which represents the mRNA expressionlevels of n genes across a number of experiments. Some of the genes inthe data have known labels according to a functional class. Fivefunctional classes were selected along with genes that belong uniquelyto these five functional classes. This yielded a dataset with 208 genes,with 79 features (experiments). Data was normalized by subtracting themean and dividing by the standard deviation for each column. This wasalso performed for the rows, and repeated for the columns. At this stagethe first three principal components were extracted. The distributionand histogram of scores is give in FIG. 2 a fork over the range of 2 to7. The same behavior is observed as seen in the mixture of fourGaussians data of FIG. 1 a. Between k=5 and k=6, there is a transitionfrom a distribution that has a large component near 1, to a widedistribution that is very similar to the distribution on the randomdata. The clustering solution that was obtained for k=5 agreed well withthe given labels, with a correlation score of 0.95.

Example 3 Uniformly Distributed Data

The results of a test run on data uniformly distributed on the unit cubeis shown in FIGS. 3 a and 3 b. The distributions are quite similar toeach other, with no change that can be interpreted as a transformationfrom a stable set of solutions to unstable solutions.

The preceding examples indicate a simple way for choosing k as the valuewhere there is a transition from a score distribution that isconcentrated near 1 to a wider distribution. This can be quantified,e.g., by an increase in the area under the cumulative distributionfunction or by an increase inS(K)=P(s>0.90).The value of 0.9 is arbitrary, but any value close to 1 would work onthe set of examples considered here.

Example 4 Isolet Letter Pronunciation

The next test was run on a portion of the ISOLET (Isolated Letter SpeechRecognition) database created by Ron Cole and Mark Fanty of theDepartment of Computer Science and Engineering, Oregon GraduateInstitute, Beaverton, Oreg. and available from the UCI (University ofCalifornia at Irvine) Repository of Machine Learning Databases. (Thisdata set was generated by having 150 subjects speak the name of eachletter of the alphabet twice, generating 52 training examples from eachspeaker.) This test provides an example of what occurs with there iscluster overlap.

One thousand points of the original dataset, representing the letters“a”, “c”, “d”, “e”, “f” were used. The columns of the data were“whitened”, then the first three principle components were extracted.FIG. 4 a is a plot of the distribution of the first two principlecomponents for each of the selected letters. FIG. 4 b provideshistograms of the correlation score for each k in the range of 2 to 7,and FIG. 4 c is an overlay of the cumulative distribution of thecorrelation score with k=2 towards the far right side of the plot(near 1) and k=10 towards the left hand side of the plot.

Table 1 provides a comparison of the number of clusters identified usingthe inventive method against the number of clusters obtained using othermethods for selecting k. These other methods are described by R.Tibshirani, G. Walther, and T. Hastie in “Estimating the number ofclusters in a dataset via the Gap statistic”, Tech. Report, Departmentof Statistics, Stanford University, 2000; also published in JRSSB 2000,where the Gap Statistic methods is shown to be superior to a number ofother methods. TABLE 1 Jain Silhouette KL Gap Subsamp 4 Gauss 4 4 9 4 4(FIG. 1) Microarray 3 2 4 6 5 (FIG. 2) Isolet 2 5 5 — 3 (FIG. 4)

In the case of clustering algorithms whose output depends on the initialcondition, e.g. k-means, a distribution of scores exists even whenconsidering a fixed sub-sample. In such cases the method produces anindication of how varied the solutions are for various values of k. Wehave observed that for a “good” value of k a similar transition occurs,but is generally “smeared”, since k-means produces widely varyingsolutions. To address this, we implemented a version of k-means similarto that presented by P. Bradley and U. Fayyad in “Refining initialpoints for k-means clustering” (in J. Shavlik, editor, Proc. of the 15thInter. Conf of Machine Learning (ICML '98), pages 91-99, San Francisco,Calif., 1998. Morgan Kaufmann.). That method produces initial conditionsthat converge to near optimal solutions, and use a fixed initialcondition for each value of k. In contrast, k-means clustering accordingto the present invention produces solutions that are highly stable withrespect to sub-sampling. This good result may be due to the globaloptimization criterion, which differs from the local bottom up approachof hierarchical clustering that appears to be less stable tosub-sampling.

Example 5 Gene Expression Data

A data set of 1600 genes with 12 time steps was utilized to illustratethe process undergone by gene expression profiles.

First, the genes were ranked in order of “quality” to pre-select asubset for further analysis. All the genes were ranked according tothree criteria: (1) saliency (the absolute difference between their minand max value; the larger the better); (2) smoothness (a coefficientassessing the smoothness of the profile was computed; the smoother thebetter); and (3) reliability (the average standard deviation of theexperiment replicated for the whole profile).

The ranks of the genes according to these three criteria were then addedto form a combined criterion according to which the genes were rankedagain. The result can be seen in FIGS. 5 a and 5 b, which show geneexpression temporal profiles. The ten best genes are depicted in FIG. 5a, while the ten worst genes are shown in FIG. 5 b according to acombined criterion of saliency, smoothness, and reliability.

The 5% top quality genes according to the above defined combinedcriterion (800 genes) were selected. A kernel clustering algorithm basedon k-means was then run on random subsets of 100 genes among these 800genes. A maximum number of clusters of ten were used, but only five didnot degenerate. The stability of the solutions was verified by runningagain kernel k-means on the resulting cluster centers. The solution wasrobust with respect to increasing the number of genes (doubling, to 1600genes), changing the subset size (to 200 genes) and the maximum numberof cluster (to 20 genes).

FIGS. 6 a and 6 b illustrate the average profiles of the clusters of theeight clustering runs and their groupings into meta-clusters 1-5. Thecluster centers in FIG. 6 a were obtained with multiple runs of k-meansusing random subsets of 100 genes in the top 800 best quality genes.FIG. 6 b shows the average cluster centers for the five clusters. Only asubset of the nine possible profiles that could occur are represented.

According to the present invention, the clustering algorithm is basedon, but is a variation of, the classical k-means algorithm. Thealgorithm operates by the following routine:

-   -   Initialize: Start with a random assignment of class labels to        the patterns.    -   Step 1: Compute cluster centers by averaging class members in        each class.    -   Step 2: Re-assign patterns to the cluster with nearest cluster        center.    -   Iterate step 1 and 2 until the assignment of patterns to classes        remains constant.

This algorithm differs from the classical k-means in that it uses aspecial metric to measure proximity in step 2. It is based on theresiduals of a fit of one profile onto another, using an affinetransformation that is a combination of translation, scaling, androtation (to first order). Such fit is remarkably simple to implementwith a couple of lines of Matlab® (The MathWorks, Inc., Natick, Mass.)and is quite fast. For example, to fit a profile (a 12 dimensionalvector x2) that goes through zero between time step 5 and 6, ontoanother vector x1, one can write:

-   -   xx2=[x2,ones(12,1),[−5;−4;−3;−2;−1;1;2;3;4;5;6;7]];    -   w=xx2\x1;    -   x2fit=xx2*w;    -   residual=mean((x1−x2fit).{circumflex over (α)}2);

FIGS. 7 a-d depict two profiles and variations on how they can be fitone onto the other or both on their average. This provides andillustration of curve fitting with affine transformations. FIG. 7 ashows the two original profiles P1 and P2. FIG. 7 b shows the profile P1fitted on P2. FIG. 7 c shows profile P2 fitted on P1. FIG. 7 d showsboth profiles fitted to their average.

After some experimentation, the following measure of dissimilarity wasadopted:

-   -   residual0+max(residual1,residual2),        where residual0 is the squared Euclidean distance between x1 and        x2, residual1 is the residual of the fit of x1 onto x2, and        residual2 is the residual of the fit of x2 onto x1. The        rationale behind this choice is that a dissimilarity simply        based on the symmetric fit to the average is, in some cases, too        optimistic—it becomes very easy with the type of affine        transformations that are being allowed to fit any curve to a        straight line (but not vice versa). Residual 0 is added to avoid        allowing too large transformations. The same desirable        properties could be achieved by solving an optimization problem        under constraints to limit the range of transformations, but        this would be computationally more expensive.

Any positive monotonic transformation of the dissimilarity does notaffect the algorithm. It should also be noted that in Step 1 of thealgorithm, a simple average of the patterns is used, as opposed to anaverage of the curves fitted to the cluster centers. After iterating,there is a significant distortion of the cluster center profiles, someof which just become straight lines.

FIGS. 8 a-h illustrate the results of one run of the algorithm on asubset of 100 gene expression profiles with a maximum cluster number of10, where expression intensity is plotted with time. The clustering isobtained from the top 800 best quality genes. FIG. 8 a shows theoriginal profiles with random class label assignments. FIGS. 8 c, 8 eand 8 g illustrate fitted profiles of cluster members at successiveiterations. FIGS. 8 b, 8 d, 8 f, and 8 h show the cluster centers atsuccessive iterations.

The present invention provides a method for model selection inclustering. Most other approaches in the prior art are based on eitherthe sum squared distances within clusters, or some combination ofwithin-cluster distances and between-cluster distances. Not only dothese prior art methods impose the notion of how a cluster should look(e.g., “compact”), most of the methods performed poorly in practice. Theinventive method, on the other hand, provides a description of thestability of a clustering method with respect to re-sampling or noise.It is believed that this stability captures the notion of “inherent”structure that is the goal in the validating clustering solutions.

Another advantage of present invention is that it is useful inidentifying the absence of structure in the data. This differs from mostclustering validation algorithms of the prior art, with the exception ofthe gap statistic, which can only determine the number of clusters ifthat number is larger than 1.

1. A method for clustering data in a dataset comprising: selecting aplurality of granularity levels k; applying a clustering algorithm to asub-sample of the dataset so that k clusters are produced; for each k,selecting a plurality of sub-samples of the dataset; selecting aplurality of pairs of sub-samples calculating a similarity between theplurality of pairs; determining a distribution of the similarity betweenthe plurality of pairs; comparing the distributions for all k; andselecting as the optimum granularity level, the k corresponding to thetightest distribution.
 2. The method of claim 1, wherein the clusteringalgorithm is A k-means algorithm.
 3. The method of claim 1, wherein theclustering algorithm is hierarchical clustering.
 4. The method of claim1, wherein the step of selecting the highest granularity level comprisesselecting a number of clusters for which there is a transition from adistribution that is peaked near 1 to a wide distribution.
 5. The methodof claim 1, further comprising, after determining the distribution,calculating a cumulative distribution function for each granularitylevel, wherein the step of selecting the highest granularity levelcomprises identifying an increase in the area under the cumulativedistribution function.
 6. The method of claim 1, wherein the datacomprises gene expression coefficients.
 7. A method for analysis of datain a dataset comprising: selecting a plurality of granularity levels ina clustering algorithm; for each granularity level: inducing aperturbation in the data; selecting a plurality of pairs ofperturbations; computing a pairwise similarity for each plurality ofpairs; determining a distribution of the pairwise similarity; selectingthe highest granularity level having the tightest distribution.
 8. Themethod of claim 7, wherein the step of inducing comprises grouping thedata into a plurality of sub-samples.
 9. The method of claim 8 whereinthe step of grouping comprises applying a k-means clustering algorithmto the data.
 10. The method of claim 8 wherein the step of groupingcomprises applying a hierarchical clustering algorithm to the data. 11.The method of claim 7, wherein the step of inducing a perturbationcomprises initialization of the k-means algorithm.
 12. The method ofclaim 7, wherein the step of inducing a perturbation comprises addingrandom noise to the data.
 13. The method of claim 7, wherein the step ofselecting the highest granularity level comprises selecting a number ofclusters for which there is a transition from a distribution that ispeaked near 1 to a wide distribution.
 14. The method of claim 7, furthercomprising, after determining the distribution, calculating a cumulativedistribution function for each granularity level, wherein the step ofselecting the highest granularity level comprises identifying anincrease in the area under the cumulative distribution function.
 15. Themethod of claim 7, wherein the data comprises gene expressioncoefficients.
 16. A method for clustering data comprising a plurality ofstructured patterns, the method comprising: (a) selecting a clusteringalgorithm based on a dissimilarity measure between pairs of structuredpatterns; (b) randomly assigning class labels to the structuredpatterns; (c) defining a plurality of cluster centers by averagingstructured patterns within each labeled class; (d) measuringdissimilarity between each structured pattern and its correspondingcluster center by measuring a residual of a fit of one structuredpattern onto another structured pattern; (e) reassigning structuredpatterns to the labeled class with the most similar cluster center; and(f) repeating steps (c) through (e) until assignment of structuredpatterns to the labeled classes remains constant.
 17. The method ofclaim 16, wherein step (d) comprises using a fit that is invariant withrespect to affine transformations.
 18. The method of claim 17, whereinthe affine transformations comprise a combination of translation,scaling and rotation.
 19. The method of claim 16, further comprisingcalculating an average structured pattern, wherein the residual is anaverage residual of the fit of each structured pattern to the averagestructured pattern.
 20. The method of claim 16, wherein step (d)comprises measuring a Euclidean distance and a residual of a fit of onestructured pattern onto another structured pattern.
 21. The method ofclaim 16, further comprising, after step (c), clustering the clustercenters.
 22. The method of claim 16, wherein the clustering algorithm isa k-means algorithm.
 23. The method of claim 16, wherein the structuredpatterns are temporal profiles.
 24. The method of claim 23, wherein thetemporal profiles are gene expression profiles.
 25. A method forclustering patterns in a dataset comprising: selecting a plurality ofgranularity levels k; selecting a clustering algorithm adapted forproducing a number of cluster centers equal to k and for assigningpatterns to clusters corresponding to the cluster centers; for eachgranularity level k: (a) inducing perturbations in the dataset togenerate a modified dataset; (b) applying the clustering algorithm tothe at least one modified dataset to produce k clusters under each ofthe perturbations; (c) creating a new dataset comprising the clustercenters identified in step (b); (d) applying the clustering algorithm tothe new dataset using the same value of k clusters; (e) determining thestability of the clusterings at each granularity level k by measuringdissimilarity between data in the new dataset and the cluster center forthe cluster into which the data was assigned; measuring fit of the datato the cluster centers for all k; and selecting as the optimumgranularity level the k corresponding to the best fit.
 26. The method ofclaim 25, wherein the perturbations comprise a combination of one ormore of sub-sampling the dataset, changing initialization of theclustering algorithm, and adding noise to the dataset.
 27. The method ofclaim 25, wherein the dataset comprises gene expression profiles.